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ABSTRACT 

Perhaps  the  most  versatile  and  efficient  method  for  inferring  the  particle  size  distribution  function 
(PSDF)  of  aerosol  clouds  from  remote  light  scattering  measurements  is  the  constrained  linear  inversion 
procedure  due  to  Philips  and  Twomey.  However,  conventional  numerical  implementations  of  this  procedure 
are  subject  to  the  following  two  problems:  (1)  an  appropriate  discrete  approximation  must  be  chosen  for 
the  PSDF  which  adequately  achieves  the  correct  balance  between  the  conflicting  requirements  of  resolution 
of  detail  in  the  PSDF  and  of  efficiency  in  the  computation  of  the  solution  and  (2)  a  proper  value  must  be 
selected  for  the  Lagrange  multiplier  (smoothing  parameter)  that  adequately  reflects  the  tradeoff  between  the 
fidelity  to  the  observed  optical  data  and  the  smoothness  of  the  solution.  Consequently,  an  adequate  recovery 
of  the  PSDF,  based  on  the  constrained  linear  inversion  procedure,  is  usually  achieved  only  after  a  certain 
amount  of  tedious  preliminary  exploratory  analysis. 

An  alternative  implementation  of  the  constrained  linear  inversion  procedure  is  presented  which  over¬ 
comes  the  problems  associated  with  the  conventional  implementation.  Firstly,  an  explicit  analytical  (continu¬ 
ous)  representation  for  the  solution  of  the  constrained  linear  inversion  procedure  is  developed  which  obviates 
the  need  to  obtain  a  discrete  approximation  for  the  PSDF.  Secondly,  an  objective  procedure,  based  on  the 
principle  of  generalized  cross-validation,  is  utilized  for  the  selection  of  the  proper  value  for  the  Lagrange 
multiplier.  Taken  together,  these  two  developments  provide  the  basis  for  an  objective,  fully  automated  im¬ 
plementation  of  the  constrained  linear  inversion  technique.  Numerical  examples  of  PSDF  inversions,  obtained 
using  the  proposed  automatic  retrieval  algorithm,  are  presented  for  synthetic  aerosol  optical  extinction  and 
scattering  data. 
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INTRODUCTION 

1.  The  determination  of  the  distribution  of  particle  sizes  in  a  polydispersed  aerosol  cloud,  based  on  a 
remote  sensing  methodology,  does  not  permit  a  direct  measurement  of  the  size  spectra.  Rather,  the  infor¬ 
mation  concerning  the  particle  size  distribution  function  (PSDF)  must  be  inferred  from  the  light  scattering 
data  (viz.,  measurements  of  the  way  the  aerosol  cloud  scatters  electromagnetic  radiation)  which  provide 
either  the  spectral  extinction  as  a  function  of  the  incident  wavelength  and/or  the  scattered  light  intensity 
as  a  function  of  the  scattering  angle.  Given  some  combination  of  these  measurements,  the  recovery  of  the 
PSDF  of  the  aerosol  cloud  is,  in  essence,  an  inverse  or  indirect  problem. 

2.  The  retrieval  of  the  aerosol  particle  size  spectra  from  a  finite  set  of  imprecise  optical  scattering  data 
is  an  inherently  ill-posed  problem  in  the  sense  that  the  solution  is  both  nonunique  and  unstable  (viz.,  the 
solution  does  not  depend  continuously  on  the  data).  In  view  of  the  fact  that  a  knowledge  of  the  PSDF 
would  aid  in  the  understanding  of  the  various  physical  and  chemical  mechanisms  and  processes  which  are 
responsible  for  the  aerosol  microstructure  and  would  permit  the  development  of  more  efficient  and  reliable 
particle  diagnostic  methods  for  the  evaluation  and  assessment  of  airborne  particulate  hazards,  a  considerable 
research  effort  has  been  directed  to  the  development  of  algorithms  for  the  reconstruction  of  the  aerosol  size 
distribution  from  optical  scattering  data. 

3.  The  earliest  attempt  at  the  specification  of  a  PSDF  retrieval  algorithm  can  be  attributed  to  Yamamoto 
and  Tanaka  [I],  who  formulated  a  numerical  inversion  algorithm  based  on  the  constrained  linear  inversion 
procedure  originally  developed  by  Philips  [2]  and  Twomey  [3]  in  the  United  States  and,  independently,  by 
Tikhonov  [4]  in  the  Soviet  Union.  The  method  of  constrained  linear  inversion  has  subsequently  been  applied 
to  the  inference  of  particle  size  distributions  utilizing  a  myriad  of  light  scattering  measurements,  including 
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optical  depth  data  obtained  as  a  function  of  wavelength  [5],  backseat tered  radiance  measured  as  a  function 
of  wavelength  [6],  a  combination  of  optical  extinction  determined  as  a  function  of  wavelength  and  scattered 
light  intensity  determined  as  a  function  of  scattering  angle  [7],  and  spectral  light  distribution  and  polarization 
distribution  obtained  as  a  function  of  wavelength  and  scattering  angle  [8]. 

4.  Without  a  doubt,  the  constrained  linear  inversion  procedure,  which  utilizes  the  Philips  and  Twomey 
second  derivative  smoothing  constraint,  is  the  most  popular  procedure  for  the  inference  of  aerosol  size  distri¬ 
butions.  However,  in  the  implementation  of  this  procedure,  it  is  necessary  to  adopt  a  discrete  representation 
for  the  PSDF  /(r)  as  well  as  a  difference  approximation  for  the  second  derivative  smoothing  functional. 
Almost  invariably,  f(r)  is  approximated  by  a  piece- wise  constant  representation  whereby  the  size  domain 
(range  in  radii  r  of  the  particles)  is  partitioned  into  a  finite  number  of  intervals  over  which  /(r)  is  assumed  to 
be  constant.  This  discretization  leads  to  the  following  problem:  there  is  no  objective  procedure  for  selecting 
the  correct  number  of  discrete  segments  to  employ  in  the  approximation  of  /(r).  Ideally,  the  number  of 
segments  must  be  large  enough  so  that  the  resulting  binning  does  not  lead  to  a  loss  of  detail  in  the  PSDF 
and  at  the  same  time  not  so  large  that  the  resulting  discretized  problem  cannot  be  solved  economically 
on  a  computer.  Moreover,  the  constrained  linear  inversion  solution  depends  critically  on  the  selection  of  a 
proper  value  for  the  Lagrange  multiplier  (smoothing  parameter)  7  that  determines  the  degree  of  smoothing 
to  impose  on  the  solution.  Again,  there  is  no  objective  method  for  the  selection  of  an  appropriate  value  for 
7.  For  the  most  put,  attempts  at  a  resolution  to  this  problem  have  been  based  on  trial-and-error  procedures 
and/or  on  empirical  rules  derived  cither  from  computational  experiments  or  from  a  sensitivity  analysis  [7,9]. 

5.  In  view  of  the  problems  associated  with  the  application  of  the  constrained  linear  inversion  procedure, 
Curry  and  Kiech  [10,11]  have  proposed  two  alternative  procedures  for  the  retrieval  of  the  PSDF  from  light 
scattering  data.  The  first  procedure,  referred  to  as  the  constrained  eigenfunction  expansion,  expresses  the 
PSDF  in  terms  of  the  eigenfunctions  of  a  covariance  operator  constructed  from  the  Mie  scattering  kernels. 
The  coefficients  in  this  expansion  are  determined  so  that  the  recovered  PSDF  is  as  close  to  some  trial  solution 
(function)  as  possible  while  maintaining  the  required  consistency  with  the  noisy  data.  The  constrained 
eigenfunction  expansion  procedure  is  actually  a  generalization  of  an  analytic  inversion  procedure  developed 
by  Capps,  Henning  and  Hess  [12].  It  is  important  to  note  that  although  the  constrained  eigenfunction 
expansion  method  does  not  require  a  discrete  approximation  for  the  PSDF,  it  still  nevertheless  requires 
the  proper  selection  of  the  value  of  the  Lagrange  multiplier  7.  Curry  and  Kiech  proposed  that  7  be  found 
by  application  of  the  Residual  Relative  Variance  method,  a  process  which  leads  to  a  somewhat  complex 
double  iterative  computational  sequence  for  the  solution.  The  second  procedure,  designated  as  the  nonlinear 
regression  method,  is  precisely  as  the  name  implies — the  PSDF  is  postulated  to  possess  a  particular  functional 
form  (e.g.,  log-normal)  and  the  set  of  parameters  associated  with  this  form  is  then  estimated  directly  from 
the  data  based  on  the  least-squares  method.  This  procedure  can  yield  very  good  results  provided  the  proper 
parametric  model  is  chosen  for  the  PSDF  and  the  correct  number  of  parameters  are  utilized.  However,  it 
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should  be  emphasized  that  the  nonlinear  regression  method  can  also  lead  to  very  poor  results  if  the  latter 
two  criteria  are  not  adequately  satisfied.  Indeed,  in  view  of  the  fact  that,  in  most  cases,  there  is  no  a  priori 
information  on  the  correct  parametric  form  to  use  for  the  PSDF,  the  nonlinear  regression  procedure  cannot 
be  recommended  for  general  use  (i.e.,  for  “blind”  inversion).  Moreover,  the  nonlinear  regression  method  is 
computationally  intensive  due  to  the  nonlinear  occurrence  of  the  parameters  in  the  assumed  parametric  (i.e., 
model)  size  distribution  function. 


6.  Despite  the  introduction  of  the  constrained  eigenfunction  expansion  and  the  nonlinear  regression 
procedures,  the  constrained  linear  inversion  method  still  remains  by  far  the  most  widely  employed  method 
for  inferring  particle  size  distributions  from  optical  data.  This  is  because  the  latter  method  only  requires 
that  the  unknown  solution  satisfy  a  weak  a  prior  constraint  (i.e.,  a  smoothness  constraint)  unlike  the  con¬ 
strained  eigenfunction  expansion  and  nonlinear  regression  methods.  Moreover,  the  constrained  linear  in¬ 
version  method  can  lead  to  very  good  results  provided  a  certain  amount  of  careful  preliminary  exploratory 
analysis  is  undertaken  to  select  the  proper  discrete  representation  for  the  PSDF  and  the  proper  value  for  the 
smoothing  parameter  7.  In  view  of  this,  it  would  be  desirable  to  develop  a  fuily  automated  implementation 
of  the  constrained  linear  inversion  procedure  which  (I)  does  not  require  a  discrete  representation  for  the 
PSDF,  thus  obviating  the  need  to  select  the  proper  number  of  discrete  segments  with  which  to  divide  the 
size  domain,  and  (2)  incorporates  an  objective  method  for  the  selection  of  the  appropriate  value  for  7. 


7.  The  major  feature  that  distinguish  this  work  from  that  of  previous  authors  on  the  implementation 
of  the  constrained  linear  inversion  procedure  is  the  emphasis  on  the  development  of  a  completely  automatic 
implementation  of  the  procedure  that  requires  little  or  no  1  ser  intervention.  Consequently,  the  purpose  of  this 
paper  is  to  show  how  to  achieve  a  completely  automatic  implementation  of  the  constrained  linear  inversion 
procedure.  To  this  end,  it  is  first  shown  that,  for  a  fixed  value  of  the  smoothing  parameter  7,  an  analytical 
representation  for  the  solution  of  the  constrained  linear  inversion  method  can  be  obtained  by  utilizing 
standard  techniques  from  the  theory  of  abstract  Hilbert  spaces.  This  representation  expresses  the  PSDF  as 
a  continuously  differentiable  function  of  the  radius  r  and,  as  such,  obviates  the  requirement  for  the  selection 
of  an  appropriate  discrete  representation  as  is  required  in  the  usual  application  of  the  constrained  linear 
inversion  method.  Consequently,  the  problems  stemming  from  the  proper  discretization  of  the  PSDF  are  no 
longer  an  issue.  Next,  it  is  shown  how  the  principle  of  generalized  cross-validation  can  be  applied  to  provide 
an  objective  procedure  for  the  selection  of  an  appropriate  value  of  the  Lagrange  multiplier  (i.e.,  smoothing 
parameter).  In  this  way,  the  second  major  problem  associated  with  the  application  of  the  constrained  linear 
inversion  procedure  can  be  resolved.  Indeed,  the  coupling  of  the  analytic  representation  of  the  PSDF  with  the 
cross-validated  selection  of  the  Lagrange  multiplier,  provides  the  basis  for  a  fully  automatic  implementation 
of  the  constrained  linear  inversion  method  with  no  need  for  a  preliminary  (exploratory)  analysis  on  the  user’s 
part.  - - 
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MATHEMATICAL  FORMULATION  OF  THE  PROBLEM 

8.  Electromagnetic  theory  dictates  that  optical  scattering  data  from  a  polydispersed  aerosol  cloud 
consisting  of  an  ensemble  of  homogeneous  spherical  particles  is  related  to  the  particle  size  distribution 
function  through  a  Fredholm  integral  equation  of  the  first  kind  [13,14] 

«(*)=  I  K(x,r)f{r)dr,  (la) 

J n 

provided  that  the  aerosol  cloud  is  optically  thin  so  that  only  single  scattering  occurs  and  that  the  particles 
comprising  the  cloud  are  randomly  separated  so  that  only  incoherent  scattering  occurs  (i.e.,  there  is  no  sys¬ 
tematic  phase  relation  between  the  radiation  scattered  by  the  individual  particles).  In  eq.  (la),  s(x)  denotes 
the  theoretical  (model)  optical  data  (e.g.,  optical  extinction,  spectral  turbidity,  scattered  light  intensity,  etc.) 
determined  as  a  function  of  the  wavelength  *  =  A  and/or  the  scattering  angle  x  ~  6.  The  kernel  function 
K(x,r),  corresponding  to  the  datum  s(x),  represents  the  scattering,  backscatter,  or  extinction  cross  sections 
for  a  spherical  particle  of  radius  r.  This  function,  which  can  be  evaluated  using  Mie  theory,  is  also  dependent 
on  the  complex  index  of  refraction  m  of  the  particle.  In  the  following  treatment,  it  is  assumed  that  m  is 
known.  Finally,  /(r)  in  eq.  (la)  denotes  the  particle  size  distribution  function  and  fi  =  [/?] ,  specifies 
the  radii  limits  R\  and  Rj  (i.e.,  size  domain)  for  /(r). 

9.  A  particular  experiment  will  yield  M  numerical  values  of  «(x)  (e  g.,  obtained  for  a  finite  and  discrete 
set  of  measurement  variables  \Xi  }i=1)  which  are  contaminated  with  noise.  Consequently,  the  observational 
model  relating  the  measured  data  D;  to  the  model  data  s(xj)  is  given  by 

A  =  s(xi)  +  e„  *  =  1,2, ,  Af,  (16) 

where  e*  is  assumed  to  be  a  zero-mean  measurement  error  with  variance  <rf.  The  inverse  problem  associated 
with  eqs.  (la)  and  (lb)  can  now  be  stated  as  follows:  given  a  finite  set  of  inaccurate  optical  scattering  data 
and  the  errors  {<7*}*^  in  their  measurement,  determine  the  PSDF  /(r). 

10.  The  most  important  property  to  note  with  regard  to  this  inverse  problem  is  that  it  is  inherently 
ill-posed.  This  translates  into  the  observation  that  a  small  random  perturbation  on  the  optical  data  s(x)  will 
induce  unacceptably  large  variations  on  the  solution  /(r).  Consequently,  a  practical  solution  to  the  problem 
can  only  be  obtained  if  the  ill-posedness  is  directly  confronted  and  addressed.  If  this  is  not  done,  the 
computation  of  the  solution  from  practical  data  (cf.  eq.  (lb))  will  invariably  lead  to  wild  oscillations  in  the 
PSDF.  These  oscillations,  which  are  merely  artifacts  in  the  solution  that  have  no  direct  physical  significance, 
are  the  manifestations  of  the  ill-posed  nature  of  the  problem.  To  formulate  a  well-posed  problem,  Philips 
[2]  and  Twomey  [3]  have  proposed  that  a  regularized  solution  (or  “quasi-solution”)  /7(r)  be  obtained  by 
minimizing  a  smoothing  functional  of  the  form 

Ry(f)  =  RLs(f)  +  7Rs(f),  (2a) 
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where 


n  _  1  y' (A -*(*«)) 


(26) 


is  the  usual  weighted  least-squares  discrepancy  (misfit)  function  with  the  «-th  weight  chosen  to  be  the  inverse 
of  the  variance  of  associated  with  the  *-th  datum  A  and 


Ks,/)  =  /»($)  * 


(2c) 


is  a  stabilizing  functional  (penalty  term)  which  measures  the  degree  of  smoothness  (regularity)  in  the  un¬ 
known  PSDF.  Note  that  Rs(f)  discriminates  against  steep  curvatures  in  /(r)  and,  in  this  sense  at  least,  the 
PSDF  obtained  by  minimizing  Rs(f)  subject  to  the  constraint  that  f(r)  reproduce  the  data  to  within  some 
expected  accuracy,  must  generate  the  “flattest”  solution  consistent  with  the  available  data. 


11.  In  eq.  (2a),  y  is  a  non-negative  regularization  parameter  (Lagrange  multiplier).  The  smoothing 
functional  is  composed  of  the  linear  combination  of  the  data  misfit  Rts(f)  and  the  roughness  measure  Rs(f), 
with  the  latter  weighted  by  the  regularization  parameter  y.  This  parameter  must  be  chosen  to  reflect  the 
degree  of  smoothing  to  be  imposed  on  the  required  solution  and,  indeed,  as  y  increases,  the  roughness  term 
in  eq.  (2a)  is  given  an  increasing  significance  (i.e.,  weight)  at  the  expense  of  the  misfit  term.  Ideally,  the 
regularization  parameter  must  be  chosen  to  reflect  some  desirable  compromise  between  the  fidelity  with 
which  the  PSDF  regenerates  the  data  as  measured  by  the  misfit  term  and  the  smoothness  desired  in  the 
PSDF  as  measured  by  the  roughness  term.  It  should  perhaps  be  noted  that  the  choice  of  y  is  equivalent 
to  the  choice  of  the  misfit  tolerance  T,  which  reflects  the  degree  of  consistency  between  the  theory  and 
observations  in  the  form  Ri,s(f)  <  T.  Of  course,  T  (or,  equivalently  y)  must  be  chosen  to  adequately  reflect 
the  magnitude  of  the  noise  perturbing  the  data. 


12.  The  problem  posed  in  eq.  (2a)  is  usually  solved  [2,3,15]  by  dividing  the  size  domain  0  =  [i?i,  R2] 
into  N  non-intersecting  equal-sized  intervals  Ar;  =  rj  —  r;_;  (j  =  1,2,...,  )  with  /(r)  assumed  to  be 

constant  over  each  interval,  viz. 


/(»•)  =  />,  rj-i  <r<rj  (j  =  1,2,...,N), 

where  ro  =  R\  and  =  R2.  The  discretization  then  proceeds  as  follows.  First,  the  integral  of  eq.  (la)  is 
expressed  as  a  summation  by  utilizing  some  appropriate  quadrature  approximation  to  give 

A  /  =  s 


where  A  is  the  M  x  N  matrix  of  quadrature  coefficients  u>;/if(xtlr,)  (i  =  1,2, . . .  ,M;  j  =  1,2, . . . ,  N),  w,  are 
the  weights  associated  with  the  quadrature  rule,  /=  .. .  ,/\)T  and  8=  (s(ati),  *(*2).  •  ••  .s(*jm))T 

Second,  the  roughness  measure  Rs{f)  is  replaced  by  the  following  second  difference  approximation  : 

«=E(/f+i-2/f+/H)2- 

j=i 
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With  these  approximations,  the  smoothing  functional  of  eq.  (2a)  can  now  be  expressed  as  the  following 
discrete  representation: 

||Wd  -  WA/ 1||  +  -r/TH/,  (3) 

where  W  is  the  M  x  M  diagonal  matrix 

W  =  diag{l/<xi,  1/<t2, 


13.  The  solution  embodied  in  eq.  (4)  requires  that  an  appropriate  discretization  be  chosen  for  the  PSDF 
/(r).  Indeed,  the  value  of  N  must  be  chosen  large  enough  so  that  there  is  no  loss  of  fine  structure  in  /(r). 
On  the  other  hand,  although  /(r)  can  be  represented  to  an  arbitrary  degree  of  accuracy  by  increasing  the 
value  of  N,  this  incurs  the  disadvantage  of  increasing  the  dimension  of  the  (WA)T  (WA)  matrix  (cf.  eq.  (4)), 
thus  making  it  more  costly  to  invert  on  a  computer.  Hence,  practical  constraints  imposed  by  the  capabilities 
of  the  computer  determine  the  upper  bound  on  the  value  of  N.  Moreover,  it  is  important  to  emphasize 
that  certain  kernel  functions  K(x,r)  such  as  the  backscatter  cross  section  are  extremely  oscillatory  ar>d 
erratic  in  form  [14],  with  the  result  that  a  large  value  of  N  would  be  required  in  order  to  achieve  a  good 
quadrature  approximation  for  «(*)  (cf.  eq.  (la)).  Under  these  circumstances,  it  is  necessary  for  N  to  be 
large  and,  consequently,  it  is  necessary  to  deal  with  the  inversion  of  large  dimensional  matrices.  In  light  of 
these  problems,  it  would  be  desirable  to  obtain  an  analytical  representation  for  the  solution  of  eq.  (2a)  and 
thus  do  away  with  the  need  to  implement  a  discrete  approximation  for  /(r).  In  so  doing,  discrete  values 
of  /(r)  will  only  be  required  for  the  quadrature  approximation  of  the  model  optical  data  and  not  for  the 
representation  of  the  PSDF  per  se.  The  methodology  for  achieving  this  particular  goal  is  presented  in  the 
next  section. 
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ANALYTICAL  REPRESENTATION  OF  THE  SOLUTION  FOR  THE  CONSTRAINED  LIN¬ 


EAR  INVERSION  PROCEDURE 

14.  An  explicit  closed-form  analytical  representation  for  the  form  of  the  PSDF,  which  minimizes  the 

smoothing  functional  R-,{{)  of  eq.  (2a),  can  be  obtained  by  application  of  certain  standard  techniques  of 
Hilbert  space  theory  [16].  from  this  perspective,  it  is  convenient  to  visualize  each  admissible  particle  size 
distribution  function  /(r)  as  an  element  of  a  real  Hilbert  space  R  consisting  of  all  real-valued  functions 
A(r)  defined  on  ft  =  (i.e.,  the  size  domain)  such  that  (1)  h1  =  dh/dr  is  absolutely  continuous  and 

(2)  h"  =  d2h/dr 2  is  an  element  of  C2(Q),  the  space  of  Lebesque  square-integrable  functions  on  Q.  It  is 
convenient  to  equip  the  Hilbert  space  R  with  the  following  inner  product: 

(u,v)  =u(R1)v(Rl)  +  u'(Rl)v'(Ri)+  [  u"(r)v"(r)dr, 

Jn 

where  u,  v  €  W.  Note  that  the  norm  of  an  element  in  R  is  then  given  by  ||u||  =  (u,  u)1/2.  A  brief  comment  is 
in  order  concerning  the  present  choice  of  R  and  the  associated  inner  product.  In  order  to  penalize  undesirable 
undulations  in  /(r),  the  constrained  linear  inversion  procedure  utilizes  the  integral  of  the  second  derivative 
of  the  PSDF  over  the  size  domain  as  the  penalty  functional  (i.e.,  stabilizing  functional).  In  view  of  this,  H 
should  be  chosen  as  the  class  of  all  functions  which  possess  a  well-defined  second  derivative  and  the  resulting 
norm  (actually  the  norm  inherited  from  the  inner  product)  should  be  chosen  so  that  it  adequately  reflects 
the  “size”  of  the  PSDF  with  regard  to  its  “roughness”.  In  this  manner,  the  minimization  in  R  using  the 
given  (i.e.,  chosen)  norm  then  results  in  the  penalization  of  the  undesirable  oscillatory  behavior  in  /(r)  in 
an  entirely  natural  fashion.  With  these  definitions,  the  linear  functional  relationship  between  the  theoretical 
datum  s(xj)  and  the  PSDF  /(r)  (cf.  eq.  (la))  can  be  viewed  now  as  a  continuous  (bounded)  linear  functional 
Li  on  R  so,  by  virtue  of  the  Riesz  Representation  Theorem  [16],  there  must  exist  some  gt  6  R  such  that 

s(xi)  =  Li(f)=  f  K(zilr)f(r)dr=(gi,f)  («  =  1,2,  — , JW).  (5) 

J  n 

The  elements  gt  are  called  the  representers  for  the  bounded  linear  functionals  Li  (i  =  1,2, ... ,  M). 

15.  To  find  an  analytic  representation  for  the  solution  to  eq.  (2a),  it  is  convenient  to  form  two  subspaces 
from  R.  First,  it  is  useful  to  form  the  finite-dimensional  subspace  spanned  by  the  representers  gi  (i  = 
1,2,...,  A/)  of  the  problem  and  define  Q  =  sp {ffi.Sa,. ..  ,?Af}  where  sp{}  denotes  the  span  of  a  set  of 
elements  in  R.  Second,  an  inspection  of  the  stabilizing  functional  Rs(f)  (cf.  eq.  (2c))  reveals  that  all 
elements  of  R  which  cause  this  functional  to  vanish  (viz.,  Rs(f)  =  0)  must  belong  to  a  two-dimensional 
subspace  of  R  spanned  by  the  two  elements 

hi(r)  =  1  and  hj(r)  =  r  —  Ri . 

Consequently,  it  is  useful  to  define  K  =  «p{  Ai ,  A2 } .  Note  that  elements  from  K.  are  ideal  particle  size 
distribution  functions  in  the  sense  that  they  are  optimally  smooth  with  respect  to  the  stabilizing  functional — 
viz.,  any  element  k  €  K,  results  in  Rs{k)  =  0  and,  consequently,  these  elements  are  not  penalized  by  the 
stabilizing  functional. 
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16.  Nov, ,  for  a  given  value  of  the  smoothing  parameter  7,  finding  the  function  /7(r)  6  which  minimizes 
the  smoothing  functional  R-,(f)  is  equivalent  to  finding  the  element  of  Ji  which  is  “closest”  to  the  subspace 
1C  of  ideal  functions  and  is,  at  the  same  time,  consistent  with  the  data  misfit  constraint  embodied  by  Rls(/) 
It  is  convenient  to  write  the  solution  /7  in  the  form 

1  M 

/t  =  +  +  <T' 

,=  1  •= 1 

where  a  is  an  element  of  Ti  which  is  orthogonal  to  both  Q  and  £,  viz.  (<r,<fr)  =  0  (i  =  1,2,. ..  ,M)  and 
(<r,hj)  —  0  (j  =  1,2).  While  the  form  for  /7  is  quite  intuitive,  it  should  be  noted  that  this  form  follows 
directly  from  the  Decomposition  Theorem  [16]  for  Hilbert  spaces.  Note  that  the  first  term  on  the  right  hand 
side  of  eq.  (6)  is  simply  a  representation  of  a  general  element  from  1C.  Similarly,  the  second  term  on  the  right 
hand  side  of  eq.  (6)  is  a  representation  of  a  general  element  of  Q.  The  element  a  then  represents  that  part  of 
/7  in  H  which  does  not  lie  in  either  1C  and/or  Q.  It  will  now  be  shown  that  a  must  vanish  if  /-,  of  eq.  (6)  is 
to  minimize  eq.  (2a).  Since  /7  is  the  “closest”  element  to  K  that  is  consistent  with  the  data  constraints,  this 
requires  that  the  length  of  the  element  p  =  a<fl*  +  a  be  minimized  subject  to  the  data  constraints  (cf. 
eq.  (6)).  However,  note  that  only  the  first  term  X^i=i  a'9'  P  can  a^ect  the  data  misfit  since  the  model  data 
s(*,)  are  determined  only  by  the  the  ji’s  (cf.  eq.  (5))  and  ( o,gi )  =  0  («  =  1,2, ... ,  M).  Consequently,  as  far 
as  the  data  misfit  constraint  is  concerned,  <r  can  be  freely  chosen.  However,  for  minimum  ||p||,  this  requires 
that  <r  =  0.  As  a  result,  the  PSDF  /7(r)  which  minimizes  Ry(f)  for  a  fixed  value  of  the  regularization 
parameter  7  must  possess  the  form 

M  2 

A  (r)  =  5Z<*<?<(r)  +  5^kA(r). 

•=1  y=i 

17.  The  undetermined  coefficients  aj  (i  =  1,2,..., M)  and  b,  ( j  =  1,2)  in  the  representation  of  eq. 
(7)  can  be  determined  as  follows.  If  eq.  (7)  is  substituted  into  eqs.  (2a),  (2b),  and  (2c),  it  is  relatively 
straightforward  to  show  that  the  smoothing  functional  Ry(f)  can  be  expressed  as  the  quadratic  form 

/£,(/)  =  M-1  (d-  Ta  -  Sb)TW3(d-  ra-Sb)  +  yaTTS,  (80) 


where  T  is  the  M  x  M  Gram  matrix  defined  as 

r  =  (r.r)!:';22;:" .  r„  =  (*.»>, 

S  is  the  M  x  2  matrix  defined  as 

s  =  ($>)Li121...j#  -  S*i  -  (fb’hy). 

and  a  and  6  are  M-  and  2-tuples,  respectively,  given  by 

a  =  (ai,a2,  ...,aM)T  and  b=(bi,b2)T. 


(86) 


(8c) 


(8d) 
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The  data  vector  d  and  the  diagonal  matrix  W  have  already  been  defined  with  reference  to  eq.  (3).  Now, 
differentiating  Ry(f)  in  eq.  (8a)  with  respect  to  the  coefficient  or  parameter  vectors  a  and  6  and  setting  the 
results  to  zero,  yields  the  following  system  of  linear  equations: 

(r  + Af7W-J)a  +  S6=  d,  (9a) 

ST  a  —  0;  (96) 


or,  equivalently,  in  matrix  form 


r  + M7W- 
ST 


o)  (!)=©' 


where  0  is  a  2  x  2  zero  matrix  and  0  is  a  2  x  1  zero  vector.  The  solution  of  eq.  (9c)  provides  the  values  of 
the  unknown  coefficients  and  {6j}^_r 

18.  Explicit  expressions  for  the  coefficient  vectors  a  and  6  can  be  derived  as  follows.  First,  note  that  the 
matrix  (T  +  MyW~2)  is  nonsingular  since  (1)  the  Gram  matrix  T  is  symmetric  and  non-negative  definite 
by  definition  (cf.  eq.  (8b)),  (2)  W  is  positive  definite,  and  (3)  7  is  positive  by  assumption  (i.e.,  recall  that  7 
is  the  positive  regularization  parameter).  Consequently,  if  eq.  (9a)  is  premultiplied  by  ST(r+  My  W~2)  \ 


one  obtains 


sTa  +  sT(r  + a/7W-2)~‘s6  =  sT(r  +  M7W-2)-'< 


which  can  be  expressed  as 


6=  Sr(r  +  M7W' 


J)  's  sT(r 


+  MyW~2y1dt 


on  using  eq.  (9b)  and  assuming  that  the  matrix  S  has  full  rank.  Now,  a  can  be  computed  from 

a=  (r  +  A/7W"2)_1(d-S6)  (106) 

by  virtue  of  eqs.  (9a)  and  (10a).  Eqs.  (10a)  and  (10b)  provide  the  coefficients  of  the  desired  solution 
to  the  constrained  linear  inversion  problem.  Indeed,  eqs.  (10a)  and  (10b)  in  conjunction  with  eq.  (7) 
provide  a  continuous  representation  for  the  PSDF  /7(r)  which  minimizes  R-f(f).  However,  to  be  able  to 
utilize  this  continuous  representation,  it  is  necessary  to  develop  expliv.it  expressions  for  the  representers  gt 
(i  =  1,2, ... ,  M).  The  solution  of  this  problem  is  the  subject  of  the  next  section. 


EXPLICIT  COMPUTATION  OF  THE  REPRESENTERS  FOR  THE  CONSTRAINED  LIN¬ 
EAR  INVERSION  PROCEDURE 

19.  Although  the  Riesz  Representation  Theorem  guarantees  that  every  bounded  linear  functional  L, 
defined  on  a  Hilbert  space  H  possesses  a  representer  y,-  €  H  such  that  L.;(f)  =  (gi,f)  for  every  /  £  W,  the 
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theorem  does  not  state  how  such  a  function  can  be  constructed  in  general.  Hence,  it  is  useful  to  approach  the 
problem  of  the  construction  of  the  representers  for  the  constrained  linear  inversion  procedure  in  the  following 
heuristic  manner.  To  this  end,  it  is  convenient  to  view  the  Hilbert  space  H  (defined  in  the  previous  section) 
as  the  direct  sum  of  two  Hilbert  spaces  Ho  and  Hi  (viz.,  H  —  Ho®Hi)  defined  as  follows: 

He  =  {/  :  df/dr  is  absolutely  continuous,  d2f  /dr2  =  oj 

and 

Hi  =  {/  :  df/dr  is  absolutely  continuous,  d2f/dr 2  £  C2(Cl),f(R\)  =  0 ,df(Ri)/dr  =  oj. 

Observe  that  Ho  is  a  two-dimensional  subspace  of  H  which  coincides  with  the  subspace  K  defined  earlier. 
The  inner  products  on  Ho  and  Hi  must  be  defined  as 

(/.*)  o  s  f(Ri)g(Ri)  +  f‘(RiW(Ri).  f,g£H0 


(f>  9)1  =  [ 
J  n 


_  f  d2f(r)  d2C( r) 


dr 2  dr2 


dr,  f,g£H  1, 


respectively,  so  as  to  be  consistent  with  the  inner  product  chosen  for  H  (cf.  previous  section  for  the  definition 
of  this  inner  product).  Then,  (f,g)  =  {Pof,Pog)o  +  {P\f,P\g)i  ( f,g  €  H),  where  P0  and  Pi  denote  the 
projection  operators  from  H  onto  Ho  and  Hi,  respectively. 

20.  As  the  first  step  in  the  construction  of  the  representers  y,-,  consider  the  problem  of  finding  functions 
k(rV  £  Ho  and  €  Hi  such  that 


/(r')  =  /) 0.  for  any  /  £  H0 


and 

f(r')  =  for  any  f  €  Hi-  (116) 

In  the  preceding  equations,  r'  is  an  arbitrary  fixed  point  in  Q.  It  should  be  noted  that  and  Jfc^  are  simply 
the  representers  for  the  point  evaluation  functionals  in  Ho  and  Hi,  respectively.  Due  to  the  simplicity  of  the 
inner  product  on  Ho  and  in  view  of  the  fact  that  Ho  is  a  two-dimensional  space  spanned  by  the  functions 
hi(r)  =  1  and  =  (r  —  ^i),  a  little  reflection  shows  that  the  function  k^V  possesses  the  form 

klV(r)=l  +  (r-Ri)(r'-Ri),  r,r'  £  D.  (12) 


That  this  form  for  k^°\r)  is  correct  can  be  readily  verified  by  substitution  into  eq.  (11a). 

21.  The  function  fc^(r)  is  more  difficult  to  determine,  but  can  be  found  as  follows.  First,  note  that  eq. 
(lib)  can  be  expressed  as 


(13a) 
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(136) 


where  G^r)  is  a  Green’s  function  determined  such  that 

g(/.r,)=o 

where  6(  )  denotes  the  Dirac  delta  function.  The  solution  of  the  preceding  boundary  value  problem  (i.e.,  eq. 
(13b))  can  be  easily  shown  to  be  given  by 

r'l-i  /  °i  if  R\  <  r  <  r';  /14s 

That  eq.  (14)  is  valid  can  be  easily  verified  by  substituting  this  expression  into  eq.  (13a)  and  integrating 
the  result  by  parts  twice.  Now,  in  light  of  eq.  (13a),  ib^(r)  is  related  to  the  Green’s  function  as 

*#\r) 


dr 2 


=  G(r',r) 


or,  equivalently,  as 


(15) 


*'.l)(r)  =  f  G(r't  r")G(r,  r")  dr". 

Jn 

Inserting  eq.  (14)  into  eq.  (15)  and  evaluating  the  result  yields 

'  (y  —  r")(r  —  r")  dr"  =  (r  —  f?1)2(r/  —  Ry)/2  —  (r  —  R\)3/6,  if  Rt<r<r'; 

<  jR'<  (16) 

/  (rl-r")(r-r")dr"={r-Rl)(r,-Rl)2/2-(r'-Rl)3/6,  if  Ri  >  r  >  i> . 

Ur, 


#V)  =  ^ 


22.  Having  determined  the  representers  for  the  point  evaluation  functionals  in  Wo  and  Ji\ ,  it  is  now  a 
simple  matter  to  compute  the  representers  for  an  arbitrary  bounded  linear  functional  in  7i.  In  particular,  one 
is  interested  in  the  representers  gx  associated  with  the  data  functionals  Li  defined  in  eq.  (5)  (s  =  1,2,...,  M). 
To  this  purpose,  first  note  that  by  virtue  ofeqs.  (11a)  and  (lib),  /(r7)  =  (k^!\Pof)o  +  {k$\P\f)i  =  (kr’ ,  /) 
with  kr>  =  Jfcl?!  +  k[V  (/  6  W).  Consequently, 

s(®i)  =  (9i<  f)  —  Li(r> ) (/(r* ))  =  Lj(r<)(fcr« ,  /)  =  (Lj(r>)irr' , /), 

<7i(r)  =  Lj(r.)(fcr.(r))  =  [  I\(xi,r')kr'{r)  dr' , 

Jn 


so 


(17a) 


where 

*,,(r)  =  *<?>(»■) +  #V)  (17*) 

and  L,(r')(  )  denotes  the  »-th  data  functional  (cf.  eq.  (5))  operating  on  what  follows  (viz.,  (•))  considered  as 
a  function  of  r'.  Now,  substitution  of  the  functional  forms  for  k^V  and  k^.V  exhibited  by  eqs.  (12)  and  (16), 
respectively,  into  eq.  (17a)  leads  to  the  following  explicit  expression  for  the  representers  gx: 

*(r)=  f  K{xi,  r')  dr'  +  (r  -  W.)  /  (r'  -  Rl)K(siy)dt' 

J  n  Jn 

+  Ur-Rl)  f  (r,-R1)2K(xiy)dr'-1-  f  (r7  -  R1)3K(xi,  r')  dr' 

J  Jrx  d  JRi 

+  \{r-Rlf  j*\r' -Ri)K{xitr')dr' -\{r-Ri)3  f  ’  K(xi,r')dr'.  (18) 


UNCLASSIFIED 


UNCLASSIFIED  /12 

23.  With  the  analytical  form  for  the  representers  gi  given  in  eq.  (18),  it  is  now  possible  to  compute  the 
analytical  solution  for  the  constrained  linear  inversion  procedure  using  eqs.  (7),  (8b),  (8c),  (10a),  and  (10b). 
To  simplify  the  computational  procedure  somewhat,  it  is  important  to  observe  that  the  first  two  terms  of 
the  expression  for  g,(r)  given  by  eq.  (18)  are  elements  from  K  and,  consequently,  can  be  incorporated  with 
the  term  Y%mi  h^j{r)  €  AC  in  the  representation  of  /7(r).  Hence,  it  is  convenient  to  write  the  representation 
for  /7(r)  as 

M  2 

A(r)  =  Eo.«(r)  +  E»i*i(r),  (19“) 

i=i  j=i 

where 

9i(r)  =  ~(r  -  Ri)  f  (r' -  JZi)  J/C(ii,  r')  dr' -  ^  T  (r'  -  Ri)3K(xit  r')dr' 

1  JRi  6  JRt 

+  \(r-Ri  )*  J  V -Ri)K(xiy)dr'  -i(r- j  ’  K^.r'Jdr'.  (196) 

Observe  that  given  by  eq.  (19b),  is  simply  j,(r)  with  the  component  in  K  removed  (viz.,  =  P\gi). 
The  coefficient  vectors  a  and  6,  which  determine  /7(r)  in  eq.  (19a),  are  still  given  by  eqs.  (10a)  and  (10b). 
However,  in  this  case,  the  elements  of  the  Gram  matrix  F  are  given  now  by 

(gi,  9i)  =  (§i>  9j)  =  jf  dr<  »,J  =  1,2,-  -,M,  (20a) 

with 

=  jf  ^ (r'  —  Ri)K(xi,r')drJ  —  (r  —  R\)  jf  *  K(zi,r')di'.  (206) 

Furthermore,  the  elements  of  the  S  matrix  are  unchanged  and  are  given  by 

(gi,  hj)  =  Li{hj)  =  [  K(*i ,  r)hj(r)dr,  i=  1,2, ... ,  M;j  =  1, 2.  (20c) 

J  n 

24.  In  summary  then,  for  a  given  value  of  the  regularization  parameter  y,  the  solution  of  the  constrained 
linear  inversion  procedure  (i.e.,  the  particle  size  distribution  function  which  minimizes  the  smoothing  func¬ 
tional  of  eq.  (2a))  possesses  the  continuous  representation 

/7(r)  =  STg(r)  +  bTh(r),  r  e  ft,  (21a) 

where  g(r)  and  h(r)  are  vector-valued  functions  on  fl  defined  as 

g(r)  —  (gi(r),gj(r),. . .  ,gM(r))T  (216) 

and 

h(r)  =  (hl(r),hJ(r)f.  (21c) 

The  coefficient  or  parameter  vectors  a  and  6  of  eq.  (21a)  are  computed  according  to  eqs.  (10a)  and  (10b) 
with  the  matrix  elements  of  T  and  S  determined  as  per  eqs.  (20a),  (20b),  and  (20c).  This  solution  for  the 
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PSDF  (cf.  eqs.  (21a),  (21b),  and  (21c))  should  be  contrasted  with  the  conventional  approximate  form  of 
the  solution  represented  by  eq.  (4).  Eq.  (4)  is  based  on  a  discrete  approximation  (i.e.,  a  piece-wise  constant 
approximation)  of  the  PSDF  which  requires  the  selection  of  an  appropriate  number  of  discrete  size  intervals 
(N)  to  use  for  the  representation  of  the  continuous  aerosol  particle  size  distribution.  On  the  other  hand, 
eqs.  (21a)-(21c)  constitute  a  continuous  representation  of  the  PSDF  which  solves  the  constrained  linear 
inversion  problem.  With  this  solution,  there  is  no  direct  need  to  utilize  a  discrete  approximation  for  f(r) 
with  its  concomitant  subjectivity  in  the  selection  of  N .  It  is  important  to  note  that  closed-form  expressions 
cannot  be  obtained  for  the  model  optical  data  s(x<)  (cf.  eq.  (la))  or  for  the  matrix  elements  I\;  of  T  (cf. 
eqs.  (20a)  and  (20b))  and  Sij  of  S  (cf.  eq.  (20c)),  so  that  some  form  of  quadrature  rule  is  required  in 
the  computation.  However,  even  though  the  quadrature  approximation  for  s(x,)  must  necessarily  involve 
the  discretization  of  the  PSDF  /(r),  it  should  be  emphasized  that  this  discretization  is  only  utilized  to 
obtain  an  approximation  for  the  model  optical  data  and  is  not  related  to  the  representation  of  f(r)  per 
se.  In  the  usual  implementation  of  the  constrained  linear  inversion  procedure,  the  representation  of  /(r) 
as  a  piecewise-constant  function  is  directly  associated  with  the  quadrature  approximation  for  s(xi).  In  the 
present  implementation  of  this  procedure,  a  continuous  and  exact  representation  is  provided  for  /(r)  and 
a  discretization  of  this  representation  is  utilized  only  for  the  quadrature  approximation  of  a(x,) — viz.,  the 
discretization  is  not  utilized  for  the  approximate  representation  of  /(r)  as  in  the  usual  implementation. 
Finally,  the  computation  of  the  PSDF  via  eq.  (21a)  involves  the  inversion  of  a  matrix  of  order  M,  whereas 
the  computation  of  the  PSDF  via  eq.  (4)  requires  the  inversion  of  a  matrix  of  order  N.  In  the  problem  of 
the  inference  of  the  PSDF  from  light  scattering  data,  N  (number  of  size  intervals)  is  almost  always  larger 
than  M  (number  of  data  points). 

SELECTION  OF  THE  REGULARIZATION  PARAMETER 


25.  Before  a  complete  solution  can  be  given  for  the  constrained  linear  inversion  procedure,  it  is  necessary 
to  address  the  problem  of  the  proper  selection  of  the  regularization  (smoothing)  parameter  y.  In  this  paper,  it 
is  proposed  that  the  generalized  cross-validation  (GCV)  method  be  utilized  for  the  selection  of  this  parameter. 
The  GCV  method  is  a  rotation-invariant  form  of  the  ordinary  cross-validation  method  introduced  by  Allen 
[17].  Indeed,  the  GCV  method  has  been  successfully  applied  to  determining  the  correct  degree  of  smoothing 
of  discrete,  noisy  data  using  spline  functions  [18,19]  and  to  the  optimal  smoothing  of  probability  density 
[20]  and  spectral  density  estimates  [21].  For  a  detailed  theoretical  treatment  of  the  properties  of  the  GCV 
estimate  of  a  smoothing  parameter,  the  reader  should  consult  Craven  and  Wahba  [19]. 


26.  With  reference  to  the  problem  treated  in  this  paper,  the  GCV  method  specifies  that  the  optimal 
estimate  y  of  the  smoothing  parameter  y  be  chosen  to  minimize  the  criterion  function 


1  "  (Dk  -  LkifP))7  , 

* 


(22a) 


Jb- 
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where  fy ^  denotes  the  minimi  zer  of  the  smoothing  functional  Ry(f)  (cf.  eq.  (2a))  with  the  Jb-th  data  point 
(i.e.,  Du)  omitted  from  the  computation.  Recall  that  L*  (fc  =  1,2, ....  M)  are  the  data  functionals  defined 
as  per  eq.  (5).  Moreover,  in  eq.  (22a),  isj^  are  weights  given  by 


«4T)  =  (i  -p**(t)) /(l  “  j^tr(P(lr)))  ' 


where  Pkk(y)  is  the  fcfc-th  entry  of  the  M  x  M  matrix  P(-y)  defined  as 


M/r) 

MA) 
Lu  (A) 


=  P(y)d, 


and  tr(-)  denotes  the  trace  operation.  It  should  be  noted  that  the  basic  logic  behind  the  selection  of  7 
to  minimize  V(y)  stems  from  the  fact  that  this  choice  for  7  minimizes,  on  the  average,  the  squared  error 
between  the  data  points  Dk  (k  =  1, 2, . . . ,  M)  and  the  predictors  Lt  (fyk^)  of  these  points.  Observe  that  the 
predictor  for  the  fc-th  data  point  is  obtained  using  the  estimate  fyk^  of  the  PSDF  which  does  not  involve  this 
particular  data  point  (i.e.,  fyk ^  is  the  solution  to  the  constrained  linear  inversion  procedure  with  the  t-th 
data  point  omitted). 

27.  The  form  of  V{y)  given  in  eqs.  (22a)  and  (22b)  is  somewhat  complex  and  is  not  suited  for  numerical 

computations.  However,  it  can  be  shown  [19]  that  ^(7)  can  be  expressed  in  the  following  form  that  is  more 
amenable  to  computations: 

(24) 

[«-llr(I-P(7))j 

It  should  be  emphasized  that  ||  •  ||f  in  eq.  (24)  denotes  the  usual  Euclidean  norm  and  not  the  norm  ||  •  || 
inherited  from  the  inner  product  on  W.  The  remaining  problem  now  is  to  determine  the  form  of  the  matrix 
P(7)  for  the  present  problem  and  to  substitute  it  into  eq.  (24).  To  this  end,  first  note  that  (cf.  eq.  (19a)) 

M  7 

MAW)  =  E«.M*W)  +  X>;MW»-)) 

.=1  i= 1 

M  ] 

=  E0<^*-&)+  ^2bA9k,hj). 

i=l  j= 1 

In  view  of  the  preceding  expression  and  the  definition  of  the  matrix  P(7)  given  in  eq.  (23),  it  follows  that 


MA) 

MA) 

MA) 


=  T5+S6=  P(7)d, 


(l-p(7))d=  MjW~2S, 
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on  using  eq.  (9a).  In  view  of  this  equation,  it  follows  that  the  numerator  of  V(7),  as  displayed  in  eq.  (24), 
can  be  expressed  as  follows: 

||W(l  -  P(7))d||2£  =  (Af7)SaTW~2a.  (25) 

28.  To  proceed  further,  it  would  be  necessary  to  substitute  the  expression  for  a  displayed  in  eq.  (10b) 
into  eq.  (25)  and  then  try  to  manipulate  the  result  into  a  computationally  useful  form.  However,  it  turns  out 
that  the  expression  of  a,  as  embodied  in  eq.  (10b),  is  not  the  most  useful  form  for  this  purpose.  As  a  result, 
it  is  sensible  to  consider  an  alternative,  but  equivalent,  expression  for  a  that  will  prove  to  be  more  amenable 
to  simplifications,  at  least  as  far  as  the  evaluation  of  eq.  (25)  is  concerned.  Ib  this  end,  first  observe  from 
eq.  (9b)  that  a  must  lie  in  the  null  (or  annihilator)  space  of  ST  so  it  is  possible  to  express  a  as 

a  =  Qv  where  STQ  =  0  and  v  £  (26) 

In  eq.  (26),  Q  is  an  M  x  (Af  —  2)  matrix  whose  columns  span  the  null  space  of  Sr .  Again,  as  in  the  derivation 
of  eq.  (10a),  it  is  assumed  that  the  rank  of  S  is  two  (i.e.,  S  is  full  rank).  Now,  if  eq.  (9a)  is  premultiplied 
by  Qr  and  combined  with  eq.  (26),  it  is  easily  shown  that 

«  =  Q  JqT  (r  +  AfyW~2)  q]  QTd.  (27a) 

Moreover,  if  eq.  (9a)  is  premultiplied  by  STW2  and  combined  with  eq.  (9b),  the  following  expression  for  b 
results: 

6=  (sTW2s)-1STW2(<r— To).  (27  b) 

It  should  be  noted  that  eqs.  (26),  (27a)  and  (27b)  constitute  alternative,  but  equivalent,  expressions  for  a 
and  b  contained  in  eqs.  (10a)  and  (10b). 

29.  It  is  important  to  point  out  that  the  matrix  Q  is  not  uniquely  determined  by  eq.  (26).  Consequently, 
it  is  possible  to  utilize  this  fact  to  choose  a  convenient  form  for  Q  that  permits  a  simplification  in  the 
evaluation  of  ||W(l  —  P(7))d  |||..  In  this  connection,  it  is  convenient  to  choose  Q  so  that  its  columns  span 
the  null  space  of  ST  (cf.  eq.  (26))  and  satisfies 

QrW“2Q  =  I.  (28a) 

With  this  choice,  insertion  of  eq.  (27a)  into  eq.  (25)  followed  by  some  simple  matrix  manipulations  leads  to 
the  following  expression: 

||W(I  -  P(7))d|&  =  (Af7)2r(A  +  M7I)~Y 

=<*’>>  <2“> 

In  eq.  (28a),  an  eigenvector-eigenvalue  decomposition  has  been  introduced  for  QrrQ,  viz. 

QTrQ  =  UAUt,  (28c) 
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where  U  is  an  ( M  —  2)  x  (M  —  2)  orthogonal  matrix  whose  columns  are  composed  of  the  eigenvectors  of 
QTrQ  and 

A  =  diagfA^Aj,  ...(Aa/.j)7,  (28d) 

is  an  (M  —  2)  x  (M  —  2)  diagonal  matrix  whose  diagonal  elements  are  composed  of  the  eigenvalues  A, 
(i  =  1,2,...,  ( M  -  2))  of  QTrQ.  Moreover,  in  eq.  (28b),  the  vector  (  is  defined  as 

6  =  (6,6,  •  •  •  ,6tf-2)T  =  UTQrd!  (28e) 


In  a  similar  manner,  it  is  possible  to  show  that  the  denominator  of  eq.  (24)  can  be  written  as  follows: 


tr(l  -  P(y))  =  (My)tt((A  +  My  I)'1) 
li-7 


=  (My) 


1 


l A,  +  My 


(29) 


Now  substitution  of  eqs.  (28b)  and  (29)  into  eq.  (24)  finally  leads  to  the  following  expression  for  K(7): 

n)  =  /(\i  +  My)f 

i/(a.  + a#t)]2 


(30) 


In  light  of  eq.  (30),  it  is  now  relatively  easy  to  obtain  a  value  of  y  that  minimizes  V(t)  by  using  a  golden 
section  search  in  one  dimension.  This  optimal  value  for  y  can  then  be  substituted  into  eqs.  (10a)  and  (10b), 
or  equivalently,  into  eqs.  (27a)  and  (27b)  for  the  computation  of  5  and  b.  Once  this  is  completed,  the  solution 
of  the  constrained  linear  inversion  procedure  can  then  be  obtained  by  evaluation  of  eq.  (21a). 


ALGORITHM  FOR  THE  AUTOMATIC  RECOVERY  OF  THE  PSDF 


30.  In  view  of  the  results  obtained  in  the  previous  sections,  the  following  algorithm  is  proposed  for  the 
automatic  reconstruction  of  the  PSDF,  given  a  finite  number  M  of  discrete,  noisy  light  scattering  data  A 
along  with  the  standard  deviations  <r<  in  their  measurement.  A  flowchart  of  the  implementation  used  to 
realize  the  automatic  recovery  of  the  PSDF  is  given  in  Fig.  1. 

31.  Step  1:  Compute  the  matrix  elements  of  T  and  S  as  per  eqs.  (20a),  (20b),  and  (20c).  Since  the 
Mie  scattering  kernels  are  complicated  functions,  these  matrix  elements  cannot  be  evaluated  analytically. 
Rather,  the  computation  of  these  matrix  elements  must  be  effected  by  utilizing  an  appropriate  numerical 
quadrature  rule. 

32.  Step  t:  Compute  the  QR  decomposition  [22]  of  WS,  viz. 

ws  =  (q,  :  <»>)(v)’ 
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where  Rj  is  a  2  x  2  upper  triangular  matrix  and  Qi  and  Q2  are  M  x  2  and  M  x  ( M  —  2)  matrices,  respectively, 
such  that  Qf  Qi  =  I,  Q[Qj  =  0,  and  QjQj  =  I.  Now  observe  that  matrix  WQj  can  be  identified  with  Q 
in  eqs.  (26)  and  (28a)  since  QjWS  =  0  (or,  equivalently.  ST(WQj)  =  0)  and  CtfQj  =  QrW-2Q  =  I. 

33.  Step  3:  Compute  the  eigenvector-eigenvalue  decomposition  of  QTrQ  (cf.  eq.  (28c)),  construct  the 
vector  £  as  per  eq.  (28e),  and  compute  the  GCV  criterion  function  1^(7)  using  eq.  (30).  Now  apply  a  golden 
section  search  [23]  in  one  dimension  to  determine  the  non-negative  value  of  the  smoothing  parameter  7  that 
minimizes  V{y). 

34.  Step  4'  Using  the  value  of  7  determined  in  Step  3,  compute  the  parameter  vectors  a  and  b  using 
either  eqs.  (27a)  and  (27b)  or  eqs.  (10a)  and  (10b).  However,  since  Q  has  already  been  computed  in  Step 
2,  it  is  perhaps  more  convenient  to  use  eqs.  (27a)  and  (27b)  for  this  computation.  Consequently,  eqs.  (27a) 
and  (27b)  are  used  in  the  computation  for  the  present  implementation. 

35.  Step  5:  Evaluate  the  PSDF  /7(r)  (i.e.,  the  solution  of  the  constrained  linear  inversion  procedure) 
using  eq.  (21a). 

36.  It  is  important  to  note  that  the  bulk  of  the  computations  in  the  preceding  algorithm  is  contained  in 
Step  1.  However,  once  an  experimental  design  has  been  specified  (i.e.,  once  the  wavelengths  and  scattering 
angles  at  which  the  scattered  light  is  to  be  measured  have  been  selected),  it  is  only  necessary  to  compute  the 
matrix  elements  of  T  and  S  once  since  these  elements  are  independent  of  the  data.  These  results  can  then 
be  saved  for  utilization  in  the  inversion  of  all  subsequent  data  sets  that  adhere  to  the  given  experimental 
arrangement.  Generally  speaking,  the  incorporation  of  the  present  algorithm  into  an  Mie  scattering-based 
particle  sizing  instrument  would  involve  the  a  priori  computation  of  the  elements  of  T  and  S  based  on 
the  experimental  configuration  of  the  instrument.  These  matrix  elements  would  then  be  stored  in  the 
onboard  computer  and  the  implementation  of  the  PSDF  inversion  algorithm  would  then  essentially  involve 
the  computation  of  Steps  2  to  5.  Moreover,  the  computation  of  the  matrix  elements  in  Step  1  and  of  /7(r) 
in  Step  5  can  be  efficiently  implemented  using  lookup  tables  for  the  relevant  Mie  kernels.  This  would  allow 
for  a  rapid  real-time  PSDF  inversion  using  only  the  capabilities  of  modern  desktop  computers. 

SOME  NUMERICAL  EXAMPLES 

37.  If  *(x)  in  eq.  (la)  is  identified  with  the  optical  extinction  of  light  at  the  wavelength  x  =  A  for  a 
homogenous,  polydispersed  aerosol  cloud  characterized  by  the  particle  size  distribution  function  f{-)  (~r, 
equivalently,  the  normalized  number  concentration  of  aerosol  particles  in  the  size  interval  r  to  r  +  dr),  then 
the  kernel  function  K(x,  r),  corresponding  to  the  Fredholm  integral  operator  that  maps  /(r)  to  s(x),  is  given 
by 

K(x,  r)  =  xr  JQext(x,  r),  (31a) 
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where  Qn t(x,  r)  is  the  extinction  efficiency  factor  evaluated  using  Mie  theory  [13,14]  for  a  particle  of  radius  r 
and  wavelength  x.  Furthermore,  if  s(x)  is  identified  with  the  light  scattering  intensity  (i.e.,  radiance)  at  the 
scattering  angle  x  S0  obtained  for  a  fixed  wavelength  Ao,  then  the  kernel  function  K(x,r)  that  transforms 
/(r)  to  *(*)  via  the  Fredholm  integral  operator  assumes  the  form 

*(,.,)=  (£)’(ali (314) 

where  <i(x,  r)  and  ii(x,  r)  are  the  Mie  intensity  parameters  [13,14]  for  scattered  light  with  perpendicular 
and  parallel  polarization,  respectively.  It  should  be  noted  that  the  extinction  efficiency  factor  as  well  as 
the  Mie  intensity  parameters  are  complicated  functions  of  r  and  x  which  are  usually  expressed  as  infinite 
sums  of  terms  involving  the  Ricatti- Bessel  functions  and  the  Legendre  polynomials.  Nevertheless,  there  now 
exist  very  efficient  and  robust  procedures  for  the  evaluation  of  the  Mie  kernels  [14,24].  Strictly  speaking, 
it  is  important  to  note  that  the  extinction  and  scattering  parameters  also  depend  on  the  complex  index  of 
refraction  m  of  the  particles.  In  the  following  examples,  the  index  of  refraction  of  the  particles  is  assumed  to 
be  known  a  priori.  In  particular,  all  simulations  considered  are  performed  assuming  that  the  particles  have 
a  real  (viz.,  the  particles  are  non-absorbing)  refractive  index  m  =  1.54. 

38.  Optical  extinction  (or  turbidity)  data  were  calculated  at  15  discrete  wavelengths  spanning  the 
interval  between  x  =  0.1  pm  and  x  =  7.5  pm  for  various  functional  forms  of  the  PSDF  /(r).  The  radii  limits 
which  determine  the  extent  of  the  size  domain  for  /(r)  are  taken  to  be  Q  =  [Ri,  fl2]  =  [0.1  pm,  5.0pm].  In 
addition,  light  scattering  intensity  data  were  generated  at  A0  =  0.65  pm  for  15  scattering  angles  between  0  and 
90  degrees.  Since  eq.  (la)  with  the  kernel  functions  of  eqs.  (31a)  and  (31b)  cannot  be  integrated  analytically 
due  to  the  complexity  of  the  Mie  kernels,  it  is  necessary  to  calculate  the  optical  data  numerically  by  utilizing 
some  quadrature  rule.  The  present  simulations  were  generated  using  Simpson’s  rule  as  the  quadrature 
formula.  These  simulated  optical  data  were  then  degraded  with  pseudorandom  normally  distributed  noise 
adjusted  to  correspond  to  a  prescribed  level  of  root-mean-square  (RMS)  error.  These  degraded  spectral 
extinction  and  scattered  radiance  data  comprised  the  components  of  the  data  vector  d  (cf.  eq.  (lb))  of 
dimension  M  —  30  that  served  as  the  input  to  the  inversion  algorithm. 

39.  Fig  2  displays  the  result  of  the  reconstruction  of  a  uniform  PSDF  given  by  /(r)  =  1.0  for  r  €  fi. 
The  figure  illustrates  that  for  the  case  of  noise-free  multispectral  extinction  and  scattered  radiance  data,  the 
recovered  PSDF  coincides  exactly  with  the  true  (actual)  PSDF.  The  regularization  (smoothing)  parameter, 
computed  using  the  generalized  cross-validation  procedure,  was  effectively  zero  (to  within  the  roundoff  errors 
of  the  computer)  for  the  case  of  noise-free  data.  Moreover,  an  examination  of  the  coefficient  (parameter) 
vectors  returned  by  the  retrieval  algorithm  revealed  that  5=0  and  b  =  (1.0,0.0)T,  a  result  which  would 
yield  /(r)  =  1.0  exactly  as  per  eq.  (21a).  The  inversion  result  of  this  particular  example  is  interesting 
because  it  illustrates  the  fact  that  if  a  solution  can  be  constructed  from  elements  in  K  that  satisfy  the 
data  constraints,  that  solution  is  preferably  chosen  by  the  algorithm.  Recall  that  the  elements  in  K  are  the 
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preferred  solution  candidates  since  these  elements  are  considered  to  be  optimally  smooth  with  respect  to 
the  stabilizing  functional  Rs()  in  the  sense  that  Rs(f)  =  0  for  all  /  6  1C.  To  extend  the  present  noise-free 
example,  the  dashed  line  in  Fig.  3  illustrates  the  recovered  PSDF  for  the  case  when  the  input  optical  data  is 
corrupted  with  3  percent  RMS  Gaussian  noise.  For  the  case  of  noisy  data,  observe  that  a  small  component 
from  Q  has  been  selected  by  the  algorithm  in  addition  to  the  major  component  from  K.  (i.e.,  the  function 
h\(r)  =  1  from  K).  Note,  in  particular,  that  although  a  small  component  from  Q  has  been  selected  in  the 
inversion,  the  regularization  parameter  selected  by  the  generalized  cross-validation  procedure  yields  not  only 
a  relatively  smooth  solution  but  also  one  that  provides  a  good  approximation  to  the  true  solution. 

40.  Fig.  3  presents  the  recovered  aerosol  particle  size  distribution  functions  for  the  case  of  a  log-normal 
distribution  possessing  a  single  mode  at  0.5  (i m.  Two  inversions  were  performed  for  this  example:  one  for 
noise-free  optical  data  and  one  for  optical  da<.<>  that  had  been  degraded  with  5  percent  RMS  Gaussian  noise. 
Again,  the  cross-validation  procedure  selected  7  «  0  (actually  7  =  1.2  x  10-9)  for  the  noise-free  inversion 
example.  However,  for  the  noise-corrupted  example,  the  cross-validation  procedure  provided  a  value  of 
3.2  x  10~*  for  7.  Observe  that  the  inverted  PSDFs  compare  favorably  with  the  true  Bize  distribution  with 
the  correct  mode  determined  at  0.5  fim.  Again,  note  that  the  values  of  7  selected  by  the  generalized  cross- 
validation  procedure  are  appropriate  for  both  the  noise-free  and  noise-corrupted  cases  since  the  resulting 
solutions  obtained  using  these  values  appear  to  be  optimally  smooth — the  selection  of  7  is  neither  too  small 
resulting  in  an  undesirable  oscillatory  behavior  in  the  solution  nor  too  large  resulting  in  the  oversmoothing  of 
the  solution  with  the  concomitant  loss  of  relevant  detail  in  the  PSDF  (e.g.,  the  single  mode  in  the  function). 
However,  it  is  important  to  point  out  that  the  recovered  solution  for  the  case  of  noise-corrupted  d„ta  exhibited 
(cf.  Fig.  3)  negative  (and  hence,  unphysical)  values  at  the  tails  of  the  log-normal  distribution.  This  can 
probably  be  attributed  to  the  the  fact  that  the  information  content  embodied  in  the  optical  data  on  the 
nature  of  the  distribution  in  the  tails  is  small.  Consequently,  features  in  the  ends  of  the  distribution  cannot 
be  accurately  recovered  using  only  the  given  noisy  optical  data.  Furthermore,  the  magnitudes  of  the  small 
oscillations  which  result  in  the  negative  values  near  the  ends  of  the  distribution  are  below  the  level  of  noise 
corrupting  the  input  data. 

41.  Fig.  4  displays  the  recovered  particle  size  distributions  for  the  case  of  a  bimodal  log-normal  distri¬ 
bution  function  with  a  primary  mode  centered  at  0.5  ftm  and  a  secondary  mode  centered  at  2.0  fan.  The 
inversions  were  performed  on  noise-free  optical  data  as  well  as  data  that  had  been  corrupted  at  the  5  percent 
RMS  noise  level.  Note  that  the  resulting  inverted  solutions  are  good  in  that  both  modes  in  the  distribution 
have  been  properly  identified.  The  two  peaks  in  the  recovered  PSDF  fall  close  to  the  true  values  and  the  areas 
under  these  peaks  approximate  the  actual  areas  quite  well.  The  values  of  7  selected  by  the  GCV  procedure 
appear  to  provide  stable  solutions  in  both  cases,  although  it  should  be  noted  that  in  the  noise-corrupted 
case,  the  solution  exhibits  a  small  undulation  at  roughly  the  position  corresponding  to  the  relative  minimum 
between  the  two  modes.  Fig.  5  shows  the  recovered  PSDFs  for  a  Junge-like  distribution  [25]  obtained  from 
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both  noise-free  optical  data  and  data  that  had  been  degraded  with  5.0  percent  RMS  noise.  Both  inversion 
solutions  provide  good  estimates  for  the  size  distribution,  although  there  are  some  small  (minor)  oscillations 
superimposed  on  the  generally  monotonically  decreasing  trend  for  the  noise-corrupted  solution.  Fig.  6  dis¬ 
plays  the  inverted  solutions  for  the  case  of  a  Rosin-Rammler  distribution  [25].  Both  noise-free  optical  data 
as  well  as  data  corrupted  by  Gaussian  noise  at  the  5  percent  RMS  level  were  employed.  Again,  note  that 
the  agreement  between  the  inverted  and  the  actual  size  distributions  is  good.  Indeed,  the  recovered  PSDFs 
in  both  cases  reproduce  the  shape  of  the  true  distribution  well.  This  fact,  coupled  with  the  observation  that 
the  inverted  results  are  relatively  smooth,  is  evidence  that  the  Lagrange  multipliers  y  have  been  properly 
selected  by  the  GCV  method  in  both  the  noise- free  and  noise-corrupted  cases. 

CONCLUSIONS 

42.  The  purpose  of  this  paper  has  been  to  develop  an  automatic  algorithm  for  the  “blind”  reconstruction 
of  particle  size  distribution  functions  of  polydispersed  aerosol  clouds,  composed  of  homogeneous  spherical 
particles  of  known  composition  (i.e.,  of  known  refractive  index),  from  discrete,  noisy  light  scattering  data.  The 
algorithm  is  based  on  an  alternative  implementation  of  the  popular  constrained  linear  inversion  procedure 
due  to  Philips  and  Twomey.  This  procedure  tries  to  minimize  artifacts  (i.e.,  oscillations)  in  the  solution  by 
utilizing  a  “roughness”  measure  (i.e.,  stabilizing  functional)  depending  on  the  second  derivative  of  the  PSDF. 
The  automatic  algorithm  arises  from  the  attempt  to  address  two  major  problems  that  plague  conventional 
implementations  of  the  constrained  linear  inversion  procedure,  namely,  the  need  to  select  an  appropriate 
discrete  approximation  for  the  PSDF  and  the  importance  of  obtaining  a  proper  value  for  the  Lagrange 
multiplier  (i.e.,  smoothing  parameter)  that  adequately  reflects  the  amount  of  smoothing  to  impose  on  the 
solution  in  relation  to  the  level  of  noise  corrupting  the  input  data. 

43.  To  address  the  first  of  these  problems,  it  was  shown  that  an  analytic  representation  for  the  solution 
of  the  constrained  linear  inversion  procedure  can  be  obtained  by  application  of  some  standard  techniques 
from  the  theory  of  Hilbert  spaces.  This  represen tation  provides  a  convenient  parametrization  for  the  PSDF 
and,  consequently,  permits  the  estimation  of  the  entire  size  distribution  rather  than  merely  its  values  at  nodes 
corresponding  to  some  discretization  of  the  particle  size  domain.  This  differs  from  the  usual  implementation 
of  the  constrained  linear  inversion  procedure  whereby  discrete  values  are  used  both  for  the  approximation  of 
the  PSDF  and  for  the  quadrature  required  to  evaluate  the  model  optical  data.  By  constrast,  in  the  present 
implementation  of  the  procedure,  the  discrete  values  of  the  PSDF  are  utilized  only  for  the  quadrature.  The 
second  problem  was  addressed  by  application  of  the  principle  of  generalized  cross-validation,  which  provides 
an  objective  procedure  for  the  selection  of  an  appropriate  value  of  the  Lagrange  multiplier.  With  an  ob¬ 
jective  method  for  the  computation  of  the  correct  value  of  the  smoothing  parameter  and  with  the  explicit 
representation  for  the  PSDF,  the  reconstruction  of  the  aerosol  size  distribution  function  is  reduced  to  the 
problem  of  estimating  a  set  of  unknown  coefficients  whose  values  then  determine  the  exact  functional  form  of 
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the  PSDF.  This  estimation  problem  can  be  efficiently  resolved  as  the  solution  of  a  system  of  linear  algebraic 
equations.  It  has  been  demonstrated  using  synthetic  data  examples  that  the  automatic  implementation  of 
the  constrained  linear  inversion  procedure  works  well  and  requires  a  minimum  of  user  intervention.  These 
examples  seem  to  indicate  that  the  GCV  procedure  provides  reliable  estimates  for  the  regularization  pa¬ 
rameter  that  adequately  reflect  the  noise  level  corrupting  the  data.  Indeed,  the  automatic  algorithm  when 
implemented  with  a  lookup  table  method  for  the  Mie  kernels,  provides  the  possibility  for  a  rapid  real-time 
PSDF  inversion  procedure  using  standard  desktop  computers. 

44.  In  this  paper,  the  GCV  procedure  was  coupled  with  an  analytic  (i.e.,  exact)  representation  for  the 
solution  of  the  constrained  linear  inversion  procedure  in  order  to  construct  an  automatic  retrieval  algorithm 
for  the  PSDF.  However,  it  is  important  to  emphasize  that  the  GCV  procedure  can  be  coupled  with  other 
analytic  representations  for  the  PSDF  in  order  to  formulate  alternative  automated  recovery  procedures. 
In  this  connection,  the  eigenfunction  expansion  technique  of  Curry  and  Kiech  [10,11]  can  be  automated 
by  interfacing  it  with  the  GCV  procedure.  It  is  interesting  to  note  that  the  selection  of  the  regularization 
parameter  using  the  GCV  procedure  is  closely  related  to  the  the  selection  using  the  Residual  Relative  Variance 
(RRV)  procedure  proposed  by  Curry  and  Kiech.  Indeed,  both  procedures  attempt  to  select  the  regularization 
parameter  so  that  the  computed  “true”  predictive  root- mean-square  error  (i.e.,  the  residual  error  between 
the  noisy  observations  and  the  model  observations  computed  using  the  true  PSDF)  is  approximately  equal  to 
the  experimented  error.  However,  the  RRV  procedure  requires  the  implementation  of  a  somewhat  awkward 
double  iterative  calculation  procedure  whereby  a  PSDF  inversion  computation  is  required  for  every  iteration 
used  in  finding  the  optimal  value  for  the  regularization  parameter.  By  contrast,  the  GCV  procedure  decouples 
the  iterative  computation  for  the  regularization  parameter  from  the  PSDF  inversion  calculation. 

45.  The  primary  drawback  in  the  application  of  the  GCV  procedure  to  the  constrained  linear  inver¬ 
sion  technique  resides  in  the  fact  that  it  is  possible  to  sometimes  obtain  negative  values  for  the  recovered 
PSDF.  These  negative  values,  which  usually  arise  from  high-frequency  oscillations  exhibited  in  the  tails  of 
a  sharply-peaked  distribution  (viz.,  a  Gibbs  phenomenon)  is  probably  the  result  of  the  small  information 
content  present  in  the  data  concerning  the  features  in  the  tails  of  the  distribution.  There  are  two  ways  to 
overcome  the  negativity  in  these  estimates:  (1)  by  increasing  the  value  of  the  regularization  parameter  to 
such  a  point  that  there  are  no  longer  negative  values  in  the  reconstructed  distribution  or  (2)  by  provid¬ 
ing  a  representation  for  the  PSDF  that  explicitly  incorporates  the  positivity  constraints.  The  first  remedy 
provides  an  indirect  method  for  addressing  the  negativity  problem  and  has  been  applied  by  numerous  re¬ 
searchers  [5,7,9].  Unfortunately,  this  method  may  result  in  an  oversmoothing  of  the  PSDF  in  order  to  satisfy 
the  non-negativity  constraints.  Furthermore,  the  GCV  procedure  is  no  longer  applicable  to  this  case.  On 
the  other  hand,  the  second  remedy  is  more  attractive  in  that  the  non-negativity  constraints  are  directly 
applied  and  the  GCV  procedure  is  still  applicable.  In  this  connection,  the  problem  needs  to  be  solved  with 
the  constraints  a  >  0  and  6  >  0  imposed.  Unfortunately,  a  closed-form  solution  cannot  be  written  for  this 
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constrained  problem  in  contrast  with  the  unconstrained  problem  and  an  appropriate  numerical  method  needs 
to  be  developed. 
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Figure  1 


Flowchart  of  the  implementation  of  the  automatic  recovery  of  the 
particle  size  distribution  function  f(r). 
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Figure  2 

Recovered  particle  size  distribution  functions  f(r)  obtained  for  both  noise-free 
optical  data  and  data  degraded  with  3  percent  RMS  Gaussian  noise.  Also 
shown  is  the  true  (actual)  size  distribution  which  is  uniform  over  the  entire 

particle  size  domain  Q. 
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Figure  3 

Recovered  particle  size  distribution  functions  obtained  from  both  noise-free 
optical  data  and  data  corrupted  with  5  percent  RMS  noise.  The  true  size 
distribution  is  a  unimodal  log-normal  distribution  with  the  mode  centered  at  0.5  pm. 
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Figure  4 

Recovered  particle  size  distribution  functions  obtained  from  both  noise -free 
optical  data  and  data  corrupted  with  5  percent  RMS  noise.  The  true  size 
distribution  is  a  bimodal  log-normal  distribution  with  modes  centered 

at  0.5  pm  and  2.0  f*m. 
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Figure  5 

Recovered  particle  size  distribution  functions  obtained  from  both  noise-free 
optical  data  and  data  corrupted  with  5  percent  RMS  Gaussian  noise.  The  true 
size  distribution  is  a  Junge-type  distribution  shown  by  the  solid  line. 
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Figure  6 

Recovered  particle  size  distribution  functions  obtained  from  both  noise -free 
optical  data  and  data  corrupted  with  5  percent  RMS  Gaussian  noise.  The  true 
size  distribution  is  a  Rosin -Rammler  distribution  shown  by  the  solid  line. 
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Perhaps  the  most  versatile  and  efficient  method  for  inferring  the  particle  size  dis¬ 
tribution  function  (PSDF)  of  aerosol  clouds  from  remote  light  scattering  measurements  is 
the  constrained  linear  inversion  procedure  due  to  Philips  and  Twomey.  However,  con¬ 
ventional  numerical  implementations  of  this  procedure  are  subject  to  the  following  two 
problems:  (1)  an  appropriate  discrete  approximation  must  be  chosen  for  the  PSDF  which 
adequately  achieves  the  correct  balance  between  the  conflicting  requirements  of  resolu¬ 
tion  of  detail  in  the  PSDF  and  of  efficiency  in  the  computation  of  the  solution  and  (2)  a 
proper  value  must  be  selected  for  the  Lagrange  multiplier  (smoothing  parameter)  that 
adequately  reflects  the  tradeoff  between  the  fidelity  to  the  observed  optical  data  and 
the  smoothness  of  the  solution.  Consequently,  an  adequate  recovery  of  the  PSDF,  based  on 
the  constrained  linear  inversion  procedure,  is  usually  achieved  only  after  a  certain 
amount  of  tedious  preliminary  exploratory  analysis. 

An  alternative  implementation  of  the  constrained  linear  inversion  procedure  is  pre¬ 
sented  which  over-comes  the  problems  associated  with  the  conventional  implementation. 
Firstly,  an  explicit  analytical  (continuous)  representation  for  the  solution  of  the 
constrained  linear  invet sion  procedure  is  developed  which  obviated  the  need  to  obtain  a 
discrete  approximation  for  the  PSDF.  Secondly,  an  objective  procedure,  based  on  the 
principle  of  generalized  cross-validation,  is  utilized  for  the  selection  of  the  proper 
value  for  the  Lagrange  multiplier.  Taken  together,  these  two  developments  provide  the 
basis  for  an  objective,  fully  automated  implementation  of  the  constrained  linear  inver¬ 
sion  technique.  Numerical  examples  of  PSDF  inversions,  obtained  using  the  proposed 
automatic  retrieval  algorithm,  are  presented  for  synthetic  aerosol  optical  extinction  and 
scattering  data. 
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